function srf = spheresrf(r, f, c, col)
% spheresrf  - create a new sphere SRF
%
% FORMAT:       srf = spheresrf([radius, factor, center, color])
%
% Input fields:
%
%       radius      radius in BV units (mm, default: 140)
%       factor      number of faces factor, 20 * (4 ^ f) (default: 6)
%       center      center of sphere (default: [128, 128, 128])
%       color       RGBA color (default: [1, 0.7, 0.3, 1]);
%
% Output fields:
%
%       srf         valid SRF object with specified properties

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
   ~isa(r, 'double') || ...
    numel(r) ~= 1 || ...
    isinf(r) || ...
    isnan(r) || ...
    r <= 0
    r = 140;
end
if nargin < 2 || ...
   ~isa(f, 'double') || ...
    numel(f) ~= 1 || ...
    isinf(f) || ...
    isnan(f) || ...
    f < 0 || ...
    f > 9
    f = 6;
else
    f = round(f);
end
if nargin < 3 || ...
   ~isa(c, 'double') || ...
    numel(c) ~= 3 || ...
    any(isinf(c) | isnan(c) | c < -512 | c > 512)
    c = [128, 128, 128];
else
    c = c(:)';
end
if nargin < 4 || ...
   ~isa(col, 'double') || ...
    numel(col) ~= 4 || ...
    any(isinf(col(:)) | isnan(col(:)) | col(:) < 0 | col(:) > 1)
    col = [1, 0.7, 0.3, 1];
else
    col = col(:)';
end

% create output
srf = bless(BVQXfile('new:srf'), 1);
numvert = 2 + 10 * (4 ^ f);
numsurf = 20 * (4 ^ f);
srf.NrOfVertices = numvert;
srf.NrOfTriangles = numsurf;
srf.MeshCenter = c;
srf.ConvexRGBA = col;

% build arrays
vert = zeros(numvert, 3);
norm = zeros(numvert, 3);
tria = zeros(numsurf, 3);
nei = {6, zeros(1, 6)};
nei = nei(ones(numvert, 1), :);

% fill first 12 vertices neighbors
gr = (1 + sqrt(5)) / 2;
vert(1:12, :) = [ ...
   - 1,   0,  gr; ...
     1,   0,  gr; ...
   - 1,   0, -gr; ...
     1,   0, -gr; ...
     0,  gr,   1; ...
     0,  gr, - 1; ...
     0, -gr,   1; ...
     0, -gr, - 1; ...
    gr,   1,   0; ...
   -gr,   1,   0; ...
    gr,  -1,   0; ...
   -gr,  -1,   0];
nei(1:12,:) = { ...
    5, [ 5,  2,  7, 12, 10]; ...
    5, [ 1,  5,  9, 11,  7]; ...
    5, [ 4,  6, 10, 12,  8]; ...
    5, [11,  9,  6,  3,  8]; ...
    5, [ 2,  1, 10,  6,  9]; ...
    5, [ 5, 10,  3,  4,  9]; ...
    5, [11,  8, 12,  1,  2]; ...
    5, [ 4,  3, 12,  7, 11]; ...
    5, [ 5,  6,  4, 11,  2]; ...
    5, [ 5,  1, 12,  3,  6]; ...
    5, [ 2,  9,  4,  8,  7]; ...
    5, [ 7,  8,  3, 10,  1]};

% scale and move
r = r / sqrt(1 + gr * gr);
vert = [r * vert(:, 1) + c(1), r * vert(:, 2) + c(2), r * vert(:, 3) + c(3)];

% put arrays into SRF
srf.VertexCoordinate = vert;
srf.VertexNormal = norm;
srf.TriangleVertex = tria;
srf.Neighbors = nei;
srf.VertexColor = zeros(numvert, 4);
